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Abstract 

We study gravitational collapse of rapidly rotating relativistic polytropes of 
the adiabatic index T = 1.5 and 2, in which the spin parameter q = J/M 2 > 1 
where J and M are total angular momentum and gravitational mass, in full 
general relativity. First, analyzing initial distributions of the mass and the 
spin parameter inside stars, we predict the final outcome after the collapse. 
Then, we perform fully general relativistic simulations on assumption of axial 
and equatorial symmetries and confirm our predictions. As a result of simu- 
lations, we find that in contrast with the previous belief, even for stars with 
q > 1, the collapse proceeds to form a seed black hole at central region, and 
the seed black hole subsequently grows as the ambient fluids accrete onto it. 
We also find that growth of angular momentum and mass of the seed black 
hole can be approximately determined from the initial profiles of the density 
and the specific angular momentum. We define an effective spin parameter at 
the central region of the stars, q c , and propose a new criterion for black hole 
formation as q c ^ 1. Plausible reasons for the discrepancy between our and 

previous results are clarified. 
04.25.Dm, 04.30.-w, 04.40.Dg 
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I. INTRODUCTION 



One of the fundamental issues in numerical general relativity is to explore the final fate 
after gravitational collapse of rotating stellar cores. If their mass is not as large as the 
maximum mass of a neutron star ~ 2M Q , a rotating neutron star will be the outcome of the 
collapse. On the other hand, for stars of sufficiently large mass, the collapse will proceed 
completely and a spacetime singularity will be formed according to the singularity theorems 
of Hawking and Penrose [1,2]. If the cosmic censorship conjecture suggested by Penrose [3] 
is correct, any singularity should be surrounded by an event horizon. Then the black hole 
uniqueness theorems of Israel [4], Carter [5], and Robinson [6] tell that a collapsed star will 
consequently settle down to a Kerr black hole. On the other hand, if this conjecture is not 
correct, a rotating stellar collapse might form a state with naked singularities. 

It is well known that in a Kerr spacetime, the singularity is covered by an event horizon 
only if a non-dimensional spin parameter defined as q = J/M 2 , where J and M are the 
angular momentum and the gravitational mass of the system, does not exceed unity [1,2]. 
Otherwise the singularity is naked. This implies that the value of q for any black hole 
cannot be larger than unity. For the realistic progenitor of black holes, however, the value 
of q may be larger than unity (see Sec. Ill A). Thus, it is interesting to explore the final fate 
after gravitational collapse of rotating stars with q > 1. Numerical relativity is the unique 
approach to resolve this problem. 

There have been several studies in full general relativity with regard to the above sub- 
ject on the assumption of axial symmetry [7-13]. A series of simulations of rotating stellar 
collapse in full general relativity was first performed by Nakamura [7] and his collaborators 
[8] using the (2+l) + l formalism developed by Maeda et al. [14]. They adopted differen- 
tially rotating massive stars that are to collapse (i.e., whose masses are much larger than 
the maximum allowed mass for a neutron star formation). An interesting finding in their 
simulations is that q is an important parameter for determining the prompt black hole for- 
mation. Their results suggest that for q > 1, no black hole is formed and the stars bounce 
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back due to centrifugal force, indicating that the cosmic censorship conjecture holds. Stark 
and Piran [9] performed simulations for the collapse of polytropes with the adiabatic index 
T = 2 and with an artificially given rigid rotation, using the Bardeen-Piran formalism [15]. 
They reconfirmed the result found by Nakamura; q ~ 1 is the demarcation between black 
hole formation and bounce. Abrahams et al. [10] studied the collapse of axisymmetric tori 
consisting of collisionless matter, and also found that black holes form only from initial con- 
figurations with q < 1. Shibata [11] investigated the effects of rotation and shock heating 
on the criteria for prompt formation of black holes, for rotating equilibrium polytropes with 
T = 2 and q < 1, and found that effects of rotation and shock heating raise the critical mass 
for black hole formation. Shibata [12] also studied formation of black holes from marginally 
stable supramassive rotating neutron stars, modeled by T = 2.5, 2.25, 5/3, and 1.5 poly- 
tropes rotating at the mass shedding limit for which q < 1, and found that the final state of 
such collapse is a Kerr black hole with no appreciable disk. Shibata and Shapiro [13] studied 
the collapse of a rigidly rotating polytrope at the mass-shedding limit with T = 4/3. The 
value of q for such a configuration is close to unity as ~ 0.96 [16]. They indicate that the 
final state is a black hole surrounded by an appreciable disk (M disk /M 0.1). All these 
results seem to suggest that q ~ 1 is the maximum value for the black hole formation, and 
the final state after the gravitational collapse of rotating stars with q < 1 is a rotating Kerr 
black hole with a small disk mass. 

However, the rotating stellar collapse in general relativity has been not sufficiently studied 
for realistic stellar cores; e.g., iron stellar cores and pair-unstable oxygen cores for which 
T < 4/3 at the onset of the collapse. Recently, Shibata [17] determined marginally stable and 
rigidly rotating stars with soft equations of state (r < 4/3) against gravitational collapse, 
which are plausible initial conditions for rotating stellar core collapse. He indicates that (i) 
even for a rigidly rotating star with q > 1, a central region, in which an approximate local 
value of q (q c ) is smaller than unity, will first collapse to form a black hole, unless q is not 
too large (q < 2.5), and that (ii) as a result of the collapse for q c < 1 and q >> 1, a massive 
disk (M disk /M = O(0. 1)) will be formed around the final black hole. 
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In this paper, we examine the validity of the conjectures indicated in [17] by axisymmetric 
simulations in full general relativity. The essence of the points in [17] is that the density 
profile for the equilibrium stars with T < 4/3 is significantly different from that for stiff 
equations of state, only for which previous works have paid attention [8,9,11,12]. This fact 
could modify the criterion of black hole formation from previous studies. As a first step 
toward a more realistic simulation, in this paper, we focus on extracting a physical essence 
for the criterion of black hole formation using simple toy models. We perform simulations for 
rotating stellar collapse with a moderately soft equations of state of T = 1.5 and with a stiff 
equation of state of T = 2 and compare the results. From the results of these simulations, 
we illustrate that the global value of q is an adequate parameter for predicting the black 
hole formation only for the stiff equations of state but inadequate for the soft equations of 
state. 

The paper is organized as follows. In Sec. II, we briefly describe our formulation, gauge 
conditions, and boundary conditions. In Sec. Ill, we first review plausible initial conditions 
for black hole formation in which the value of q can be larger than unity. Then, we describe 
initial conditions adopted in this paper and predict the final outcomes of the stellar collapse 
for our initial models following [17]. Sec. IV presents numerical results, emphasizing that the 
predictions made in Sec. Ill are correct. Sec. V is devoted to a summary. Throughout this 
paper, we adopt the geometrical units G — c — 1, where G and c denote the gravitational 
constant and the speed of light, respectively. We use the Cartesian coordinates x k = ( x, y, z ) 
as the spatial coordinates, with r = ^/x 2 + y 2 + z 2 . 

II. SUMMARY OF FORMULATION 

We perform fully general relativistic simulations for rotating stellar collapse in axial 
symmetry using the same formulation as in [18], to which the reader may refer for details 
and basic equations. The fundamental variables for the hydrodynamics are: 

p: rest mass density, 
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e: specific internal energy, 



P: pressure, 




(1) 



where subscripts i,j,k,--- denote x, y, and z, and p the spacetime components. As the 
variables to be evolved in the numerical simulations, we define a weighted density p* = 
pcra*e 6fli and a weighted four-velocity Ui = hui = (1 + e + P/p)ui where h denotes a specific 
enthalpy. From these variables, the total baryon rest mass and angular momentum of system, 
which are conserved quantities in axisymmetric spacetime, can be defined as 



General relativistic hydrodynamic equations are solved using a so-called high-resolution 
shock-capturing scheme [19,18] on the y = plane with the cylindrical coordinates (x, z) (in 
the Cartesian coordinates with y — 0). 

To model initial conditions we adopt the polytropic equations of state 



where n is the polytropic index and K polytropic constant. Then physical units enter 
the problem only through the polytropic constant K, which can be chosen arbitrarily or 
else completely scaled out of the problem. Thus, in the following we display only the 
dimensionless quantities which are defined as 



M* = M*K~ n/2 , M = MK~ n/2 , R = RK~ n/2 , J = JK~ n , p = pK n , Q = QK n , (5) 



where M, R, and Q denote an ADM mass, a radius, and an angular velocity. Hereafter, we 
adopt the units of K — 1 so that we omit the bar. 

On the other hand, during the time evolution, we use the so-called T-law equations of 
state of the form 




(2) 



(3) 



P = Kp 



(4) 
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P = (r - l)pe, (6) 

where the adiabatic index T is set as 1 + 1/n. In the absence of shocks, no heat is generated 
and the collapse proceeds in an adiabatic manner, preserving the polytropic form of the 
equations of state. 

We neglect effects of viscosity and magnetic fields. The timescale of dissipation and 
angular momentum transport due to these effects are much longer than the dynamical 
timescale, unless the magnitude of viscosity or magnetic fields is extremely large [20]. Thus 
neglecting them is an appropriate assumption. 

The fundamental variables for geometry are: 

a: lapse function, 
/3 k : shift vector, 

7ij-: metric in 3D spatial hypersurface, 
7= e 12 ^ = det( 7 «), 

K; L f extrinsic curvature. (7) 

As in the series of our papers, we evolve 7^, 0, Aij = e~^(Kij — ^ijK k k ), and trace of 
the extrinsic curvature K k k together with three auxiliary functions Fj = 5^ k dj%k with an 
unconstrained free evolution code as in [21-25,11,18]. 

The Einstein equations are solved in the Cartesian coordinates. To impose axisymmetric 
boundary conditions, the Cartoon method is used [26]: Assuming the reflection symmetry 
with respect to the equatorial plane, simulations are performed using a fixed uniform grid 
with the grid size N x 3 x N in (x, y, z) which covers a computational domain as < x < L, 
< z < L, and —A < y < A. Here, N and L are constants and A = L/N. In the 
Cartoon method, the axisymmetric boundary conditions are imposed at y — ±A. Details 
are described in [18] to which the readers may refer. 

As the slicing condition we impose an "approximate" maximal slicing condition in which 
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K k k ks is required [21]. As the spatial gauge, we adopt a dynamical gauge condition [27]. 
In the present work, the equation for the shift vector is written as 



where At denotes a time step in numerical computation [12]. Note that in this gauge 
condition, f3 l obeys a hyperbolic-type equation for sufficiently small value of At because 
the right-hand side of the evolution equation for Fj contains a vector Laplacian term [25]. 
The outstanding merit of this gauge condition is that we can save computational time 
significantly. It has already been found that stable simulations for rotating stellar collapse 
and merger of binary neutron stars are feasible in this gauge [24,12]. 

An outgoing wave boundary condition is imposed for hij(= %j — 5^), A^, and Fj at the 
outer boundaries of the computational domain. The condition adopted is the same as that 
described in [25]. On the other hand, for and K k k , outer boundary conditions are imposed 
as r = const and K k = 0, respectively. 

A black hole may be formed as a result of collapse. We determine the location of it 
using an apparent horizon finder developed in [28]. As the system approaches a stationary 
state, the apparent horizon will approach the event horizon. In a dynamical spacetime we 
compute the apparent horizon mass Mah which is defined as [29] 



where A denotes area of an apparent horizon. 

During the numerical simulations, conservation of mass and angular momentum and a 
violation of the Hamiltonian constraint, H crror , are monitored as code checks. Here H erroi is 
evaluated with a weighted average by as 



d t (3 k = ^ kl (F l + Atd t F l ), 



(8) 




(9) 




(10) 



where V is defined as 
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(11) 



In Eq. (11), ip = e 4 ^, w = au f , and A and R are the Laplacian and the Ricci scalar with 
respect to the conformal metric 7^. Numerical results for several test calculations, including 
stability and collapse of spherical and rotating neutron stars, have been described in [18]. 



III. INITIAL CONDITIONS AND PREDICTIONS 
A. Candidates for black hole progenitor 

In this subsection we review two realistic candidates for the progenitor of black holes 
and recall that the value of q for them may be larger than unity. 

The first candidate is a degenerate iron core in excess of the Chandrasekhar mass [30,31]. 
It is well known that stars whose initial mass are larger than ~ 8M evolve to form a core 
mainly composed of iron group elements. Because of iron being the most stable nuclei and 
no energy generation by nuclear burning, the core contracts gradually. Accordingly the 
central temperature, T c , and central density, p c , rise to be T c > 10 9 K and p c > 10 9 g/cm 3 . 
Then for the stars of initial mass larger than ~ 1OM , the photo-dissociation of iron to 
lighter elements occurs first at such hot and dense environments. The entropy of electrons 
is reduced, as density increases, by the loss of thermal energy due to the photo-dissociation. 
As a result, the adiabatic index T decreases below 4/3 and the core is destabilized [30,31]. 
If mass of the core is much larger than the maximum neutron star mass, a black hole will 
be formed after the collapse. 

The second candidate is a star of very low metallicity or a member of the so-called 
first stars. Recent numerical simulations of the collapse of primordial molecular clouds 
have suggested that the first generation of stars contains many massive members of mass 
M > IOOMq [32-34]. Such massive stars may evolve to form large oxygen core after the 
helium burning phase [35]. If the oxygen core is sufficiently massive, it will encounter 
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electron-positron pair creation. As a result, thermal energy that might have gone into 
raising the temperature and providing more pressure support is diverted to it, bringing its 
adiabatic index T below 4/3 to destabilize the core. Following the pair creation, the oxygen- 
to-silicon and silicon-to-nickel transformation will supply internal energy to the core, which 
could force back the collapse to explode. However, if mass of the core is sufficiently large, 
the contraction proceeds and then nickel-to-alpha phase transition will set in, nullifying all 
decelerating effect of the nuclear burnings. Eventually the core will completely collapse to 
form a black hole [35,36]. 

The gravitational collapse sets in when stability of the progenitors changes and, hence, a 
marginally stable star against gravitational collapse may be regarded as an initial condition 
for the collapse [17]. Therefore, for both candidates described above, initial conditions for 
black hole formation are marginally stable rotating stars with adiabatic constant T < 4/3. 
Note that for a configuration with T < 4/3, rotation can stabilize it against gravitational 
collapse. 

Indeed, using polytropic equations of state, a stability criterion against gravitational 
collapse can be written, including rotational effect and post-Newtonian correction as [30,17] 

M 

3r - 4 - 2(3r - h)(3 - fc(r)— < 0, for instability, (12) 

R 

where (3 is the ratio of rotational kinetic energy, T, to gravitational potential energy, W, as 

" S PT (13) 

A; is a parameter which is ~ 6.75 in the case T = 4/3 and R is a radius of spherical polytrope. 
For the progenitor of black holes described above, the compactness M/R is lager than 10 3 
and, hence, the general relativistic correction may be neglected. To analyze the effect of 
rotation, we restrict our attention to the rigidly rotating case for simplicity in the following. 
Based on a numerical analysis in Newtonian gravity, the maximum value of (3 for rigidly 
rotating stars, which is achieved when the velocity at the equatorial surface is equal to the 
Keplerian velocity (i.e., at the mass-shedding limit), with polytropic equations of state is 
approximately written as [17] 
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0.00902 + 0.124 r 
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) 



(14) 



max 
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Combining Eqs. (12) and (14), it is found that stability criterion becomes as follows: 



Therefore, the rotation stabilizes the star with 1.328 < T < 4/3. This is a crucial result. 
In a rapidly and rigidly rotating equilibrium state, the spin parameter is related to the 
compactness parameter M/R as 



On the other hand, the value of R/M for a marginally stable star becomes larger for the 
smaller value of T [17], accordingly increasing the value of q. Recall that the value of q is 
close to unity for a maximally rotating and marginally stable star with T = 4/3 [16]. Taking 
all these into account, it is likely that the value of q exceeds unity for a rapidly rotating and 
marginally stable star with T between 1.328 and 4/3. Actually, it has been shown that this 
is the case for marginally stable and rigidly rotating polytropes with 1.328 < T < 1.332 [17]. 
Thus, the value of q for the rapidly rotating progenitor of black hole formation in nature is 
likely to be larger than unity. 



We are most interested in the final fate after gravitational collapse of rotating marginally 
stable stars with 1.328 < T < 1.332 and q > 1. Such simulation is computationally chal- 
lenging because the dynamical scale changes by a factor of ~ 10 4 during the collapse and 
also because it is necessary to take into account a realistic equation of state for getting sci- 
entific results. Fortunately, the physics that we want to understand (i.e., criterion for black 
hole formation) can be extracted even from a simulation with appropriate "toy" models. 
Here, we adopt simplified initial conditions following Stark and Piran [9]: We first give a 
marginally stable spherical polytrope, and then, add an angular momentum artificially as 



r < 1.328, for instability. 



(15) 




(16) 



B. Method of preparing initial conditions 
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well as reduce the pressure to induce the collapse. We choose two values of T as 1.5, with 
which an insight into the collapse of stars with soft equation of state may be gained, and 
T = 2, the value that Stark and Piran adopted. 

In more detail, we prepare initial conditions in the following procedure. First, we give 
a spherical star, marginally stable against gravitational collapse, using the polytropic equa- 
tions of state with T = 1.5 or 2.0. Second, we reduce the pressure by an arbitrarily chosen 
fraction fp of its equilibrium pressure. Third, an angular momentum is artificially added 
according to 



where a and denote the lapse function and the conformal factor of the spherical polytrope, 
and Q is given by 



Here, w = ^/x 2 + y 2 and R is a parameter which controls the degree of differential rotation. 
For R — > oo, the rotation approaches the rigid rotation (in the following we refer the case 
with R — > oo as the rigid rotation case). Forth, the Hamiltonian and momentum constraints 
are reimposed by solving the constraint equations, and then the time evolution is set out. 

In Tables I and II, we list characteristic quantities for our initial conditions with T = 1.5 
and 2.0 in the present work. Rs = 6.45 is a coordinate radius of a marginally stable spherical 
polytrope with T = 1.5 in the isotropic coordinates and all the quantities are scaled to be 
non-dimensional using the relations (5). 



Using the same method as described in [17], we predict the final outcome of the collapse 
for our initial conditions with q > 1. We pay particular attention to the black hole formation 
assuming that (i) the collapse proceeds in an axisymmetric manner, (ii) the viscous angular 




(17) 




(18) 



C. Predicting the final outcome 
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FIG. 1. (a)The distribution of q(j) and (b) jisco as a function of m(j)/M(= m*(j)/M*) for 
models with f^o = 0.065 and for T = 1.5. The solid and dashed curves denote the results for rigid 
rotation and differential rotation with Rq = 1.5Rs, respectively. The filled diamonds denote values 
of m(j)/M(= m^fi/M*) and q{j) at the maximum value of jisco- 

momentum transport during the collapse is negligible, and (iii) the pressure or heating 
effects never halt the collapse. Because of the assumption that the viscous effect is negligible 
during the collapse, the specific angular momentum j of each fluid element is conserved in 
an axisymmetric system. Here, j is denned as 

j = hu v . (19) 

Then we define a rest-mass distribution m*(j) as a function of j, which is the integrated 
baryon rest mass of fluid elements with the specific angular momentum less than j: 
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FIG. 2. The same as Fig. 1 but for Qq = 0.090. The solid, dashed, and dotted curves denote the 
results for rigid rotation, differential rotation with Rq = R$, and with Rq = Rg/1.5, respectively. 



m*(j) = 2ir p*r 2 drd(cos9). 

Jj'<3 

Similarly, a specific angular momentum distribution is defined according to 

J(j) = 2tt f p*j'r 2 drd{cos6). 
Ji'<j 



(20) 



(21) 



For the following analysis, it is better to have a quasi-local gravitational mass since in 
contrast with the case of soft equations of state studied in [17], the rest mass and the 
gravitational mass are different and their difference is not negligible for models adopted in 
this paper (see Table I). However, the gravitational mass cannot be locally defined in general 
relativity. Thus, here, assuming the ratio of the quasi-local gravitational mass to the rest 
mass is uniform inside a star, we define a "gravitational mass distribution" as 
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FIG. 3. The same as Fig. 2 but for model of Vt = 0.100. 



(22) 



Note that m(j) is equal to M for a maximum value of j (hereafter j max ) and that for T ~ 4/3, 
m(j) is approximately identical with m*(j). 

From these distribution functions, we define spin parameter distributions as 



7.!./) -^ ^d q{j) - J{j) 



(23) 



These may be approximately regarded as the non-dimensional spin parameters of fluid ele- 
ments with the specific angular momentum less than j. Although it is not clear which of two 
spin parameters is better for the analysis, we adopt q(j) because of the following reasons: 
(i) since q(j) is equal to the global quantity q for j = j max , q(j) would be the better quantity 
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FIG. 4. The same as Fig. 2(a) but for model of = 0.115. 

for a large value of j; (ii) M* is always larger than M for all the models. (This is likely to 
be the case for all the stable stars.) This implies that q(j) > q*(j). As a result, for q(j) < 1, 
q*(j) is always smaller than unity which is the key value for black hole formation. 

In Figs. 1(a), 2(a), 3(a), and 4, q(j) as a function of m*(j)/M* are displayed for the 
parameters listed in Table I for T = 1.5. Figure 5 also shows the same relation for param- 
eters listed in Table II for T = 2.0. These figures indicate that the value of q(j) at the 
center of stars (hereafter denoted as q c and g* jC ) can be much smaller than unity even if the 
global values of q and g*(j max ) (hereafter denoted simply as q*) are larger than unity. An 
outstanding difference between the results for two values of T is the ratio of q/q c for the 
rigidly rotating case, which is ~ 2 for T = 1.5 while ~ 1.25 for T = 2. This results from the 
fact that for softer equations of state, the star has a more centrally-condensed structure. 

As collapse proceeds, an inner region of the stars collapses faster. This property will be 
more outstanding for the softer equations of state, with which stars have a more centrally- 
condensed structure. Taking this and the distribution of q(j) into account, we conjecture that 
an inner region of q(j) < 1 may form a seed black hole first and the black hole subsequently 
grows as the ambient fluids accrete onto it. 

Assuming that a seed black hole is formed during the collapse for models of q c < 1, 
we predict the subsequent evolution in the following manner. Here, we focus only on the 
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FIG. 5. The same as Fig. 4 but for V = 2.0. The solid, dashed, and dotted curves denote the 
results for VLq = 0.68, 0.67, and 0.65 respectively. All the models are rigidly rotating cases. 

rapidly rotating case with a large value of q > 1 but with q c < 1. Since the centrifugal force 
of the rapidly rotating progenitors is large, first, the collapse is likely to proceed along the 
rotational axis to be a disk structure. In such a disk, fluid elements of a cylindrical radius 
are likely to have the approximately identical value of the specific angular momentum j. As 
a result, the fluid elements of the same value of j will collapse in a simultaneous manner 
even if they are initially at different locations. Now, let us consider the innermost stable 
circular orbit (ISCO) around the growing black hole at the center. If the value of j of a 
fluid element is smaller than that of the seed black hole at the ISCO, jisco, the element 
will fall into the seed black hole eventually. In fact, there is a possibility that some fluids 
can be captured even for j > jisco if if is 011 a noncircular orbit. Ignoring such trajectories 
yields the minimum amount of fluids that will fall into the black hole. The value of jisco 
will change as ambient fluids accrete onto the black hole. If jisco increases as a result 
of the accretion, the more ambient fluids will fall into the black hole. This suggests that 
the evolution of mass m*(j) and angular momentum J(j) inside the seed black hole will 
approximately proceed according to the initial distribution of m*(j) and J(j) because all 
viscous effects such as angular momentum transfer are assumed to be negligible. On the 
other hand, if jisco decreases during the accretion, no more fluid will fall into the black hole, 
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and as a result, the dynamical growth of the black hole will terminate. 

To estimate the value of jisco and to predict the growth path of the seed black hole, we 
assume that the spacetime metric can be instantaneously approximated by that of a Kerr 
spacetime of mass m(j) and spin q(j). On these approximations, we can compute jisco of 
a seed black hole as [30,37,38], 



Here, m(j) and q(j) are not strictly the gravitational mass and the spin parameter. We 
suspect that they may have a systematic error of magnitude ^ |1 — M/M*| for the mass and 
< 2|1 — M/M*\ for the spin parameter. This implies that the estimation of jisco by Eq. 
(24) also includes a systematic error. The magnitude of the error for jisco ma y be < 10%. 

In Figs. 1(b), 2(b), and 3(b), we show jisco [?™(j), <?(j)] as a function of m^(j)/M^ for V = 
1.5. These figures show that for models in which the degree of differential rotation is not too 
high, jisco has a maximum (hereafter denoted as jisco:max)- Thus, we predict that the seed 
black hole will grow until j reaches jiscomax- In Table III, we show jisco:max, ra*(jisco:max), 
and ^(jiscoimax) for T = 1.5. Note that for models of Q = 0.115 in which q c > 1, we 
predict that a seed black hole is not formed. Now, let us assume that m*(jisco:max)/M* 
and <?C/isco:max) m &y be approximately regarded as the mass fraction M B n/M and spin 
parameter ^bh of the final black hole (which is in a quasi- stationary state) as 



Jisco 




(24) 



where 



risco = m(j) [3 + Z 2 - {(3 - Z 1 )(3 + Z 1 + 2Z 2 )} 1/2 ' , 
Z 1 = l+[l- q(j) 2 } V3 [{1 + g(j)} 1/3 + {1 - qU)V /3 \ , 

z 2 = MJ) +zl ■ 



M BH ~ m( jisco :max 1 5 



(25) 



qBH ~ ^(jlSCOanax)- 



(26) 
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From this approximation we immediately predict that an appreciably massive disk 
(Mftg-k/M = 0.1-0.8) will be formed after stellar collapse. Note that M B h and g B H com- 
puted above might include a systematic error because of the reason that we approximate a 
spacetime composed of a black hole and a massive disk simply as a Kerr spacetime. However, 
Jisco is likely to be determined from the local quantities near the black hole and, hence, we 
consider that our treatment may be approximately correct. 

It should be addressed that for a fixed value of q c , the corresponding value of ?(jisco:max) 
is almost independent of the differential rotation parameter R , i.e., the value of q. This 
suggests that the spin parameter of the final state of a black hole may be determined by the 
value of q c and independent of q as far as q(j) is an increasing function of j. 

In summary, we have predicted the following facts in this section: (I) inner region of 
a star with q c < 1 will collapse first to form a seed black hole even if the global value q 
exceeds unity significantly; (II) the formed seed black hole will grow as the ambient fluids 
subsequently accrete onto it, (III) evolution of the relation between the rest mass m*(j) and 
angular momentum J(j) enclosed inside a growing black hole will agree approximately with 
the initial relation between m*(j) and J (J); (IV) whether black hole is formed or not will 
be determined by the value of q c independent of the value of q; (V) the final outcome of a 
dynamical collapse for a star with q c < 1 is a black hole surrounded by a disk. 

These predictions made form the initial conditions are quite reasonable, but to confirm 
them, it is obviously necessary to perform fully general relativistic simulations. The next 
section presents the results of numerical simulation and demonstrates that the predictions 
are correct. 
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(a) f^o (b) r/p~o 

FIG. 6. Evolution of (a) the central density p* iC and (b) the central value of the lapse function 
a c for T = 1.5. In both figures, solid, dashed, and dotted curves denote the differentially rotation 
models of R s = Rq/1.5 with fi = 0.090, 0.100, and 0.115 (q c = 0.87, 0.96, and 1.08), respectively. 

IV. RESULTS OF NUMERICAL SIMULATION 

A. Results for T = 1.5 

We perform simulations for various initial models with T = 1.5 listed in Table I, varying 
grid spacing as Ax = 0.018, 0.015, 0.01125, 0.009, and 0.006 (in units of G = c = K = 1). 
In the simulations, the uniform grid is adopted and the outer boundaries along the x and 
z-axes are located at L ~ 10. Numerical simulation is performed on FACOM VPP5000 in 
the data processing center of the National Astronomical Observatory of Japan, and personal 
computers with Pentium 4 processors, each of which has 2 Gbytes memory and a 3.0 GHz 
clock. 

In. Fig. 6(a) and 6(b), evolution of the central value of (p* )C ) and the central value of 
the lapse function a c for T = 1.5 and for the differentially rotation cases with R s = Ro/1.5 
and tt = 0.090, 0.100, and 0.115. We find that for models with q c < 1 {tt < 0.100), the 
collapse proceeds almost monotonically to form an apparent horizon irrespective of values 
of q and the differential rotation parameter R . On the other hand, for all the models with 
f2o — 0.115 for which q c > 1, we find that a black hole is not formed and that the stars 
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FIG. 7. Snapshots of the density contour curves of p* and of the velocity field of (v x ,v z ) at 
selected time slices for the rigid rotation model with Qq = 0.065, The contour curves are drawn for 
p* = p a x lO" ^', with j = 0, 1, 2, ■ ■ ■ , 10 where p a = 0.011, 1, and 10 at tp 1 /^ = 0, 1.45, and 1.95. 
The thick solid circle with the radius ~ 0.2M of the last panel denotes the apparent horizon. 

experience bounces and oscillation. This results from the fact that the centrifugal force near 
the axis is too strong to form a seed black hole. 

In Figs. 7 and 8, we plot snapshots of the density contour curves of p* and of the velocity 
field of (v x ,v z ) at selected time slices for the rigid rotation model with Q = 0.065 and for 
the differential rotation model with Q = 0.115 and R = R s /1.5. As we predicted in the 
last section, the collapse first proceeds along the rotational axis (z-axis) to be a disk-like 
structure. Then, the disk collapses to the center to form a black hole for q c < 1. Irrespective 
of the models with q c < 1, the collapse and the black hole formation proceed in essentially 
the same manner. For q c > 1 [see Fig. 8], the centrifugal force is strong enough to prevent 
a black hole formation. In this case, an oscillating disk is the outcome after the collapse. 
Such disk may be unstable against a nonaxisymmetric instability. 

Even if the value of q is much larger than unity, a black hole is formed for q c < 1, as 
predicted in Sec. III. To confirm this conclusion strictly, we performed convergence tests 
varying the grid resolution for a wide range. In Figs. 9(a) and 9(b), we show H error and the 



20 




X/M X/M X / M 



FIG. 8. The same as Figs. 7 but for the differential rotation model with £1$ = 0.115 and 
Ro = Rs/1-5. The contour curves are drawn for = p a x lO" ' 5 - 7 , with j = 0, 1, 2, ■ ■ ■ , 10 where 
p a = 1 for all the time steps. 

central density p* iC with different grid spacing (Ax = 0.009, 0.015, and 0.018) for f2 = 0.065 
and R = 1.5Rs- These figures illustrate that second-order convergence is achieved. 

To accurately determine the threshold values of q or q c , it is also crucial that the rest- 
mass distribution m*(J) as a function of specific angular momentum is conserved accurately 
at least until the first formation of an apparent horizon. In Fig. 10(a) and 10(b), we compare 
the rest-mass distribution at t — and at the first formation of apparent horizon for the 
rigid rotation case with Q = 0.065 and 0.090. These figures demonstrate that the rest-mass 
distribution as a function of j is conserved well, implying that a spurious numerical transfer 
of the angular momentum is negligible. 

In Fig. 11(a) and 11(b), we show time evolution of the apparent horizon mass for rigid 
rotation models with f2 = 0.065 and 0.090. These figures indicate that the process of the 
black hole formation can be divided into two phases. One is a phase in which a seed black 
hole is formed at the central region. The other is a phase in which the seed black hole grows 
as the ambient fluids fall into it. We define the mass of the seed apparent horizon MAH, S eed 
from the location of the break of the curves for Mah(£)- From this identification, it is found 
that MAH,seed is much smaller than the total mass of the system M and M B h denned in Eq. 
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FIG. 9. Time evolution of (a) violation of the Hamiltonian constraint H erTOT and (b) the central 
value of p* (p* : c) for Qq = 0.065 and i?o = l.5Rs- In panel (a), the solid, dashed, and dotted curves 
denote the results of Ax = 0.009, 0.015, and 0.018. In panel (b), the filled diamonds and open 
triangles denote the results of Ax = 0.015 and 0.018. These figures indicate that second-order 
convergence is achieved in our simulations. 

(25). As a reasonable result, it is also found that MAH,seed is smaller for the larger value of 
Q and Rq (i.e., for the more rapidly and rigidly rotating cases). 

For a Kerr black hole of mass Mbh and a spin parameter Qbh> the irreducible mass of 
the event horizon M irr is denned as [2] 



If the final state of a black hole can be approximated as a Kerr black hole even though it 
is surrounded by an appreciable disk, Mah should asymptotically approach M irr [MBH, <?bh]- 
Assuming that M B h and ^bh ma y be evaluated from the approximate relations (25) and 
(26), M irr /M is w 0.69 and 0.34 for the rigidly rotating initial data with Q = 0.065 and 
0.090, respectively (see dotted lines in Fig. 11). We note that for the differentially rotating 
cases, the value is slightly larger, but it is not significantly different for the same value of Q . 
Thus, Mah should approach these values. Actually, the value appears to approach them. 
Unfortunately, the computations crash in ~ 20M after the formation of apparent horizons, 
due to the grid stretching around the black hole horizons, and thus, we cannot confirm this 




(27) 
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(a) j/M (b) j/M 

FIG. 10. Rest-mass distribution as a function of specific angular momentum j (a) for the rigid 
rotation models with S7o = 0.065 and (b) for f2o = 0.090. In all figures, the solid curves denote the 
initial rest-mass distributions and the filled diamonds denote the rest-mass distributions when an 
apparent horizon is first formed. The grid resolution is Ax = 0.01125 for (a) and Ax = 0.006 for 

(b) . 

prediction strictly. To carry out a simulation beyond this time, the so-called black hole 
excision techniques [39] are necessary. 

It should be noticed that MAH, SC cd is always smaller than M irr by a factor of ~ 2. This 
implies that the black hole significantly grows during the accretion phase. For the realistic 
progenitor of a black hole formation that have a centrally-condensed density profile, the 
accretion phase is likely to be much longer than the collapse timescale (denned as Pcfi 2 
where p Cj o is the central density at t — 0). At the accretion phase, the system is composed of 
a black hole and a surrounding massive disk accreting onto the central black hole. Therefore, 
this accretion phase might be associated with the gamma ray bursts [40]. 

For all the models with fi = 0.100 in which q c rs 0.91-0.96, we find that a black hole is 
formed. However, it is difficult to obtain a convergent result with the present computational 
setting because of the restriction of the grid resolution. First, recall that the grid spacing 
required to follow a black hole formation is Ax < 0.1M sced)A H (not Ax < 0.1M), with which 
a black hole horizon is covered by ~ 10 grid points. For the model with f2 = 0.100, we 
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FIG. 11. Evolution of the apparent horizon mass Mah for the rigid rotation models with (a) 
S7o = 0.065 and (b) f2rj = 0.090. In both figures the horizontal lines denote the irreducible mass 
M irr with the rigid rotation models, « 0.67M for £l = 0.065 and « 0.34M for fi = 0.090. 

found that mass of the seed black hole is M see d,AH ~ 0.1M. This implies that for this model, 
the required grid spacing is Ax < 0.1M see d,AH ~ 0.01M 0.005. This value is almost the 
same as the finest grid spacing we adopt, indicating it difficult to obtain convergent results. 

From the above results, we have confirmed the conjectures (I), (II), and (IV) suggested 
in Sec. Ill C: We have found that the quasi-local value q c ~ 1 is an approximate threshold 
for a seed black hole formation and that after its formation, the black hole grows due to 
accretion. 

To confirm the conjecture (III), we compute the evolution of the total baryon rest mass 
and the total angular momentum enclosed inside an apparent horizon (denoted as to^ahO') 
and Jah(j) respectively) at each time step and compare their relation with the initial one 
between m*(j) and J(j). In Fig. 12(a), we display evolution of tti^ahU) an d Jah(J) (the 
diamonds) together with the initial relation (the solid curve) for the rigid rotation model 
with Q = 0.065. It is found that an evolutionary track of the black hole is approximately 
determined by the initial distribution of the mass and the angular momentum, until a 
numerical error is accumulated significantly. To be more convinced of this result, we also 
perform an additional simulation for a moderately rotating star (rigid rotation with fl = 
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(a) m,(j)/M, (b) m,(j)/M* 

FIG. 12. Evolution of the baryon rest mass m* ; AH(i) and the angular momentum Jah(j) 
enclosed inside apparent horizons (filled diamonds) for the rigid rotation models with (a) = 0.065 
and (b) f^o = 0.050. The solid curve denotes the initial distribution relation of m*(j) and J(j). 
The cross denote the location of w.*(iiscO:max) and J(jiscO:max), which are the predicted final 
values of the mass and the angular momentum of a black hole. 

0.050). In Fig. 12(b), we show the evolution track of m m (j) and J(j) together with the 
initial distribution for this case. This figure shows that the initial distribution certainly 
determines an approximate evolutionary track of the rest mass and the angular momentum 
enclosed inside the apparent horizon. These results confirm the conjecture (III) of the 
previous section. This conclusion is quite natural in particular in the present case because 
we initially deplete the pressure by a significant factor to quickly form a disk-like structure, 
for which fluid elements of the same value of a cylindrical radius have almost the same value 
of j. 

To this time, we have confirmed the conjectures (I) - (IV). Now, let us consider the 
conjecture (V). To confirm it, it is necessary to continue simulations for a long time after 
formation of apparent horizons. Unfortunately, computations crash in ~ 20M after the 
formation of apparent horizons, due to the grid stretching around the black hole horizons. 
To carry out a simulation beyond this time for confirmation of the conjecture (III), the 
so-called black hole excision techniques [39] are absolutely necessary. 
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(a) tJp^o (b) ?Vp^o 

FIG. 13. Evolution of (a) the central density p* c and (b) the central value of the lapse function 
a c for T = 2. In both figures, solid, dashed, and dotted curves denote the models of = 0.68, 
0.67, and 0.65 respectively. 

B. Results for T = 2.0 

For comparison with the T = 1.5 case, we also study rotating collapse adopting a stiff 
equation of state with T = 2. The simulations are performed with N = 1000 and Ax = 0.005 
(in units of G = c = K = 1) for the initial conditions listed in Table II. 

From the simulations, we find that a black hole is formed even from the initial data sets 
with q > 1. In Figs. 13(a) and 13(b), we show time evolution of the central density p* c and 
the central value of the lapse function a c for models listed in Table II. These figures indicate 
that the rigidly rotating stars with fi < 0.67 and q < 1.16 collapse to form a black hole. 
On the other hand, for rigidly rotating stars with Q > 0.68 and q > 1.17, a black hole is 
not formed. Therefore, the threshold value of q is between these two values. It should be 
addressed that the present critical value q ~ 1.2 agrees approximately with that found by 
Stark and Piran [9]. For the critical model, the value of q c is ~ 0.92, close to unity. Since 
there is no reason to believe that a particularly large value q ~ 1.2 should be the threshold 
for the black hole formation, we propose the quasi-local value q c ~ 0.92 ~ 1 as the threshold. 

The plausible reason that the critical value of q c is somewhat smaller than unity may be 
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as follows. First note that since configuration of stars with stiffer equations of state is less 
centrally-condensed and more uniform, the collapse proceeds in a more coherent manner; 
i.e., the fluid elements collapse to a black hole rather simultaneously. Accordingly, whether 
a black hole is formed or not after the collapse is likely to be not completely determined 
by local properties of the stars. Indeed, we find that ratio of the mass of the seed apparent 
horizon to the total mass of the system for T = 2.0 is much larger than that for T = 1.5 
case; for T = 2.0, MAR>eedM^ ~ 0.35 at the critical value of q C) in contrast with that for 
T = 1.5, MAH,seed/^f ~ 0.1. This indicates that for T = 2, not q c but a value of q(j) < 1 of 
a moderately large value of j (denoted as j±) may determine a black hole formation. Since 
q(j) is an increasing function of m*(j), 1 ~ q(ji) > q c - Therefore, the threshold value of q c 
may be smaller than unity for stiff equations of state. 

Note that for a "realistic" progenitor of black hole formation reviewed in Sec. Ill A, the 
associated equations of state are soft (T < 4/3) ones. Since mass distribution of such stars 
is significantly centrally-condensed, a demarcation between complete collapse and a bounce 
will be determined by properties of the central region. Consequently the critical value of 
q c may be almost unity. Furthermore, since the ratio M*/M is almost unity for this case, 
g^c ps 1 may be also used as the threshold for black hole formation. 

Before closing this section, let us explain the reason that q ~ 1 has been believed as a 
threshold for black hole formation in the previous works. Stark and Piran used a stiff (r = 2) 
equation of state [9] as we do. As described above, such stars have rather uniform density 
distribution and the collapse proceeds in a coherent manner. Furthermore, the distribution 
of q(j) expressed as a function of m*(j) inside stars with stiff equations of state is rather 
flat. Indeed, for T = 2.0 and rigidly rotating case, the ratio of q to q c is ~ 1.25 in contrast 
with the T = 1.5 cases for which the ratio is ~ 2 (compare Fig. 5 and Figs. 3(a) for rigid 
rotation models). These two facts explain the reason that Stark and Piran concluded that 
q ps 1 is an approximate threshold for black hole formation. 

Nakamura and his collaborators performed simulation for highly differentially rotating 
stars [7,8]. The higher the degree of differential rotation becomes, the flatter the distribution 

27 



of q(j) is, as indicated in Fig. 4. As a result, the difference between q and q c becomes 
small for their initial data. This is the reason that the value q — 1 was regarded as the 
demarcation between black hole formation and a bounce. Abrahams et al. [10] performed 
simulations for toroidal star clusters. Having a toroidal configuration implies that they are 
highly differentially rotating. Thus, the value of q c is likely as large as q, and hence, the 
global parameter q may be used as the threshold. 

V. SUMMARY 

We have reported new results about black hole formation in collapse of rapidly rotating 
stars with q > 1 and with moderately soft (r = 1.5) and stiff (r = 2.0) equations of 
state, analyzing initial conditions and performing axisymmetric simulations in full general 
relativity. 

The initial conditions are given, following Stark and Piran [9]; we first give marginally 
stable spherical polytropes with r = 1.5 or T = 2.0 and then artificially add an angular 
momentum. Adopting the same method described in [17], we have predicted that (I) inner 
region in which q c < 1 will collapse first to form a seed black hole even if the global value of 
q is much larger than unity; (II) the formed seed black hole will subsequently grow as the 
ambient fluids accrete onto it; (III) evolution of the relation between the rest mass m*{j) 
and angular momentum J (J) enclosed inside a growing black hole will agree approximately 
with the initial relation between m*(j) and J(j); (IV) whether black hole is formed or not 
can be determined by q c irrespective of q; (V) the final outcome of dynamical collapse will 
be a black hole surrounded by an appreciable disk. 

To confirm these conjectures, we performed fully general relativistic hydrodynamic sim- 
ulations using a high-resolution shock capturing scheme with the T-law equations of state. 
As a result of numerical simulations, we confirmed the conjectures (I), (II), (III), and (IV). 
From these results we conclude that the previous criterion for black hole formation that no 
black hole is formed for q > 1 is no longer valid in particular for soft equations of state. As 
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a new criterion for no black hole formation, we propose the condition q c ^ 1. 

Associated with conjecture (II), we have found that the process of black hole formation 
is divided into two phases: the first phase is that a seed black hole is formed and the other 
is the accretion phase that a part of the ambient fluids are swallowed into the formed seed 
black hole. At the time of the accretion phase, the system is composed of a black hole and 
a surrounding massive disk accreting onto the central black hole. In the collapse of massive 
iron cores, such accretion phase may be associated with the gamma ray bursts [40]. 

In this paper, we have focused only on the criterion of black hole formation analyzing 
simple toy models. To obtain a scientific result that can be compared with observational data 
such as gravitational waves, however, a simulation with a realistic initial condition and a 
realistic equation of state is necessary. As we indicate in this paper, a rapidly rotating black 
hole will be formed after rapidly rotating massive stellar core collapse and pair-unstable 
collapse, if the progenitors are not extremely rapidly rotating at the onset of collapse as 
q c < 1. Formation of a rapidly rotating black hole with q ~ 1 and accretion of a large mass 
onto such black hole of a large value of q are likely to be strong burst sources of gravitational 
waves for laser interferometric detectors. Thus, performing realistic numerical simulations of 
black hole formation is an important subject for predicting the gravitational waveforms. We 
plan to attack these computations in a fully general relativistic manner extending previous 
works [41,42]. 

For q c > 1, a black hole will not be formed promptly. In such case, the collapse leads 
to formation of a self-gravitating disk or torus. They will be subsequently unstable against 
nonaxisymmetric deformation. After the nonaxisymmetric instabilities turn on, angular 
momentum will be transported from the inner region to the outer region, decreasing the 
value of q(j) around the inner region below unity. As a result, a seed black hole may be 
formed. To follow these processes, it is necessary to perform a numerical simulation without 
assuming the axial symmetry. During the collapse, the typical length scale may change by a 
factor of 10 4 from R/M ~ 10 4 to 1. To follow the collapse by numerical simulation, very large 
computational resources will be necessary and, thus, the simulation for this phenomenon will 
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be one of the computational challenges in the field of numerical relativity. 
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Pc 


^0 


Ro/Rs 


M 


M* 


J 


q 


q* 


Qc 


fp 


0.005 


0.050 


oo (rigid) 


0.510 


0.530 


0.227 


0.87 


0.81 


0.50 


0.35 


0.005 


0.065 


oo (rigid) 


0.512 


0.535 


0.300 


1.14 


1.05 


0.64 


0.25 


0.005 


0.065 


1.5 


0.501 


0.534 


0.282 


1.08 


0.99 


0.64 


0.25 


0.005 


0.090 


oo (rigid) 


0.529 


0.546 


0.444 


1.59 


1.49 


0.84 


0.1 


0.005 


0.090 


1.0 


0.519 


0.541 


0.377 


1.40 


1.29 


0.86 


0.1 


0.005 


0.090 


2/3 


0.511 


0.537 


0.319 


1.22 


1.11 


0.87 


0.1 


0.005 


0.100 


oo (rigid) 


0.538 


0.552 


0.512 


1.77 


1.68 


0.91 


0.01 


0.005 


0.100 


1.0 


0.524 


0.545 


0.427 


1.55 


1.44 


0.94 


0.01 


0.005 


0.100 


2/3 


0.514 


0.540 


0.358 


1.36 


1.23 


0.96 


0.01 


0.005 


0.115 


oo (rigid) 


0.562 


0.563 


0.640 


2.03 


2.02 


1.01 


0.01 


0.005 


0.115 


1.0 


0.540 


0.553 


0.515 


1.77 


1.68 


1.05 


0.01 


0.005 


0.115 


2/3 


0.526 


0.546 


0.319 


1.54 


1.42 


1.08 


0.01 



TABLE I. Central density p c , angular velocity of rotation axis Qq, differential rotation pa- 
rameter Rq (Rs is the coordinate radius of the spherical polytrope in the isotropic coordinates), 
gravitational mass M, baryon rest mass M#, angular momentum J, the non-dimensional parameter 
(spin parameter) q = J/M 2 , q* = J/M 2 , and q c of the initial data for models with T = 1.5 in units 
of G = c = K = 1. fp denotes a pressure reduction parameter, "rigid" implies the rigid rotation 
model. 



Pc 


n 




Rq 


M 


M* 


J 


q 


q* 


q c 


fp 


0.318 


0.65 


oo 


(rigid) 


0.179 


0.197 


0.0364 


1.13 


0.93 


0.90 


0.01 


0.318 


0.67 


oo 


(rigid) 


0.182 


0.199 


0.0383 


1.16 


0.97 


0.91 


0.01 


0.318 


0.68 


oo 


(rigid) 


0.183 


0.199 


0.0393 


1.17 


0.99 


0.92 


0.01 



TABLE II. The same as I but for models with T = 2.0. 
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n 


Rq 


JISCO:max/M 


W-*(jlSCO:max)/M* 


(KjlSCOmax) 


0.050 


oo (rigid) 


2.29 


0.953 


0.792 


0.065 


oo (rigid) 


1.70 


0.815 


0.904 


0.065 


Rs/1.5 


1.77 


0.845 


0.902 


0.090 


oo (rigid) 


0.734 


0.435 


0.979 


0.090 


Rs 


0.753 


0.454 


0.982 


0.090 


1.5R S 


0.836 


0.511 


0.984 


0.100 


oo (rigid) 


0.399 


0.256 


0.992 


0.100 


Rs 


0.330 


0.222 


0.995 


0.100 


1.5R S 


0.287 


0.201 


0.997 


TABLE III. 


. The maximum 


values of jisco, 


and corresponding values of 


?™*(.?ISCo)/Mk and 



q(jisco) for r = 1.5 initial models of Q = 0.050, 0.065, 0.090, and 0.100. 
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